library(ggplot2)

mh_cauchy <- function(n=50000, start, scale=1, proposal_sd=1){
  
  samples <- numeric(n)
  samples[1] <- start
  
  for(xi in 2:n){
    current <- samples[xi - 1]
    
    # Generate a proposal from a Cauchy distribution
    proposal_sample <- rnorm(1, mean=current, sd=proposal_sd)
    
    proposal_density <- dcauchy(proposal_sample, location=0, scale=scale)
    current_density <- dcauchy(current, location=0, scale=scale)
    
    ration <- proposal_density / current_density
    
    # Accept or reject the proposal
    if(runif(1) < ration){
      samples[xi] <- proposal_sample
    } else {
      samples[xi] <- current
    }
  }
  
  return(samples)
}

mh_samples <- mh_cauchy(n=50000, start=0, scale=1, proposal_sd=1)

# Plotting the results
ggplot(data.frame(Samples = mh_samples), aes(x = Samples)) +
  geom_histogram(aes(y = ..density..), bins = 100, fill = "lightblue", color = "black") +
  stat_function(fun = dcauchy, args = list(location = 0, scale = 1), color = "red", size = 1) +
  labs(title = "Metropolis-Hastings Sampling from Cauchy Distribution",
       x = "Sample Values",
       y = "Density") +
  theme_minimal()


my_func <- function(x) {
  if(x < -5 || x > 5) return(0)
  3 * x^2 + 5
}

my_func(3)

my_func_metropolis <- function(n=50000, start=0, proposal_sd=1){
  samples <- numeric(n)
  samples[1] <- start
  
  for(xi in 2:n){
    current <- samples[xi - 1]
    
    # Generate a proposal from a normal distribution
    proposal_sample <- rnorm(1, mean=current, sd=proposal_sd)
    
    # Calculate the acceptance ratio
    ratio <- min(1, my_func(proposal_sample) / my_func(current))
    
    # Accept or reject the proposal
    if(runif(1) < ratio){
      samples[xi] <- proposal_sample
    } else {
      samples[xi] <- current
    }
  }
  
  return(samples)
}

mh_samples_func <- my_func_metropolis(n=500000, start=0, proposal_sd=5)
# df_func <- data.frame(Samples = mh_samples_func)

# summary(df_func)

# 画图
hist(mh_samples_func, breaks = 100, probability = TRUE, col = "lightblue",
     main = expression("Metropolis Sampling from f(x) = 3x² + 5", xlab = "x"))

# 叠加真实密度曲线（标准化）
curve((3 * x^2 + 5) / integrate(function(x) 3*x^2 + 5, -5, 5)$value,
      from = -5, to = 5, add = TRUE, col = "red", lwd = 2)
legend("topright", legend = "True Density (scaled)", col = "red", lwd = 2)

# Plotting the results for my_func
# ggplot(df_func, aes(x = Samples)) +
#   geom_histogram(aes(y = ..density..), bins = 100, fill = "lightgreen", color = "black") +
#   stat_function(fun = my_func, color = "blue", size = 1) +
#   labs(title = "Metropolis-Hastings Sampling from Custom Function",
#        x = "Sample Values",
#        y = "Density") +
#   theme_minimal()


